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Abstract 



We study the interactions of a possibly dense and/or quantum degenerate 
gas with driving light. Both the atoms and the electromagnetic fields are 



° . 

', represented by quantum fields throughout the analysis. We introduce a field 



theory version of Markov and Born approximations for the interactions of 
light with matter, and devise a procedure whereby certain types of products 



. 

of atom and light fields may be put to a desired, essentially normal, order. In 

o 

the limit of low light intensity we find a hierarchy of equations of motion for 

^ , correlation functions that contain one excited-atom field and one, two, three, 

etc., ground state atom fields. It is conjectured that the entire linear hierarchy 

may be solved by solving numerically the classical equations for the coupled 

system of electromagnetic fields and charged harmonic oscillators. We discuss 

the emergence of resonant dipole-dipole interactions and collective linewidths, 

and delineate the limits of validity of the column density approach in terms of 

non-cooperative atoms by presenting a mathematical example in which this 

approach is exact. 
42.50.Vk,03.75.Fi,05.30.Jp 
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I. INTRODUCTION 



Indications of Bose-Einstein condensation (BEC) in trapped alkali metal vapors have 
been reported recently At this point all direct probing of such condensates has been 

carried out optically. Correspondingly, in anticipation of BEC and the role of light in the 
experiments, optical response of degenerate atomic gases has been the subject of active 
theoretical research already for quite some time PHT2"||. Aside from the interest in BEC, in 
an evaporatively cooled gas of alkali atoms one may for the first time have a homogeneously 
broadened, weakly interacting system at such a high density that there are many atoms in 
a cubic wavelength, pA 3 ^> 1. This kind of a sample would in its own right serve to further 
our understanding of the interactions of light with matter. 

Nonetheless, in spite of all the theoretical work, there still are quite basic unsettled issues 
in the theory of the optical properties of dense and/or degenerate gases. Under the condition 
pA 3 ^> 1 the atoms no longer respond to the electromagnetic fields individually, but their 
properties are modified by the presence of nearby atoms. For instance, the atoms exhibit 
collective linewidths and line shifts. Inasmuch as it comes to the near-resonant response 
in the regime pA 3 ^> 1, all treatments of the optical properties of degenerate gases known 
to the present authors (including ours) contain uncontrolled approximations that bear on 
linewidths and line shifts. As a result, the regions of validity and the relations between 
different treatments tend to be somewhat ill-defined. 

A rigorous study of atom-field interactions valid regardless of atom density, atom statis- 
tics, optical detuning, and so forth, is clearly called for. The paper of Morice et al. || is a 
step in this direction. They start from a full quantum mechanical Hamiltonian, including 
quantized light and internal degrees of freedom and cm. motion of the atoms. However, at 
an early stage these authors go over to a classical treatment of the cm. motion of the atoms. 
They then derive equations of motion for a few correlation functions involving polarization 
and atom density, and solve the optical response including all photon exchange between any 
pair of atoms. 
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The program carried out in Sec. [Q| of the present paper is similar to the agenda of 
Ref. P[. We start in Sec. |11 A| from our field theory version of the Hamiltonian as in Ref. [JTTJ 



amended with the atom-atom contact interaction [|T3|,|9f that derives from the Power-Zienau 
approach. The point of departure from Ref. |J is that we describe the atoms with quantum 
fields throughout. The mathematical techniques introduced in the process are analogous 
to the time honored tools in quantum optics: field theory version of the Born-Markov 
approximation (Sec. [II B|) , and procedures to move noncommuting operators to a certain 



advantageous order (Sec. [II C|) . In the present paper we complete the derivation by assuming 



the limit of low intensity for driving light. The end result in Sec. [II D| is a hierarchy of 
equations of motion for correlation functions that involve atomic polarization at one point 
and densities at 0, 1, . . . , points in space. The lowest two equations coincide with those 
given in Ref. |§. 

In the present paper the emphasis is on the structure of the theory. To gain more insights, 
we examine a few simple special cases in Sec. [Til]. We demonstrate the exquisite subtlety of 
the propagation of radiation through an atomic sample by presenting one particular set of 
assumptions that yields the standard column density results of (optically) non-cooperating 
atoms (Sec. [ill A| ), and by reviewing the density expansion of Morice et al. (Sec. |U1 rj| ). 
In the case of only two atoms, what we call cooperative linewidth and line shift emerge as 



manifestation of the dipole-dipole interaction. This is discussed in Sec. [Ill Q . Here we also 
point out that in our limit of low intensity of the driving light, the collective linewidth and 
line shift could perfectly well have been derived from classical electrodynamics of classical 
atoms (charged harmonic oscillators). In this paper we do not attempt to derive any new 
results from the hierarchy of correlation functions, but the connection to classical physics 
points to a possible future method for exact solution of the hierarchy: classical simulations 
of a system of classical atoms. A few comments to this effect are made in Sec. |U1 D[ 

Concluding remark about the possible solutions and extensions of our hierarchy are made 
in the final Sec. [TV[ Certain mathematical details concerning the divergence of the dipolar 
field and a summary of dipole matrix elements are deferred to the Appendix. 



II. MOTION OF ATOM FIELDS 



A. Basic dynamics of the fields 



We begin by recapping, reformulating, and extending the salient results of ||10|| . The 
main items of this section are the coupled evolution equations for the light and matter fields 
exemplified by Eqs. (f|) and (0). Overall, we emphasize the similarities of the theory to the 
classical electrodynamics of polarizable media. 

1. Hamiltonian 

For better or worse, in this paper we regard atoms as point dipoles. A mathematically 
rigorous treatment produces a delta function term in the field of a dipole at the position 
of the dipole, which results in a contact interaction between dipoles. For mathematical 
consistency, this time around we therefore also keep the contact interaction generated in the 
Power-Zienau transformation from the p • A to the d • E gauge . This interaction was 



ignored as presumably inconsequential in the limit of large detuning considered in Ref. |10 
but for an arbitrary detuning it may become an issue. 

The atoms have two internal energy levels, which we label g for "ground" and e for 
"excited" . We allow for angular momentum degeneracy of the energy levels, so the complete 
specification of the internal state of an atom am includes the level label a = e or g and the 
z component of angular momentum m. We assume dipole coupling of each atom to light. 

In first quantization, we add to the Hamiltonian of [|TIJ the contact interaction, the 
polarization energy 



lfp = i-5>-di (l) 
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Here d« and r« are the dipole operator and center-of-mass position operator for the z th atom. 
There are also divergent self-energies with i = j, but we ignore these as we do not attempt 
a quantitative calculation of the Lamb shift. Equation (|I|) displays a standard two-body 



interaction, which is immediately converted to second quantization. As before, the many- 
atom system is described by Heisenberg picture quantum fields ip am (rt), which obey the 
proper commutator relations. While much of our development applies to fermions as well, 
in this paper we consider only bosons explicitly. In terms of the atom fields, the additional 
polarization energy is the integral of the Hamiltonian density 

T^P — — ^2 [d m2 Mi • d A f 2mi tpl m2 4>lM 2 ' l l J eM 1 4 , gm 1 
eo m\m 2 
M\M 2 

+ - dM 2 m 2 ' djvfimi ' l PeM 2 ' l PeM 1 ' , Pgm 2 ' > Pgm 1 + ~ d m2 M 2 ' d mi Mi 

The notation d m M stands for the matrix element (gm\d\eM) of the dipole operator of one 
atom, d. We denote the energy level implicitly in such a way that a label of a Zeeman state 
with a lower-case m refers to ground state, upper-case M to the excited state. 

However, to simplify the notation as far as possible, we are going to adopt yet another 
convention that is in force unless we explicitly state otherwise. We do not write the magnetic 
quantum numbers explicitly. For instance, we write ip e2 in lieu of ip e M 2 - Also, we write the 
matrix elements d m u as d ge . Finally, if the same level index appears twice in a product, 
sum over the magnetic substates of the level is implied. With these conventions, we write 
Eq. (^0 anew as 

Hp = —[d g2ei ■ d e2gi i)l 2 ipl 2 i) ei ip gi + \ d e292 ■ d eigi ^ViV^Vv 

+ 2 d S2e 2 ■ d 9i ei VjaV&V'eaV'ei] • (3) 



2. Electromagnetic fields 

Unlike in [|1(J], and similarly to [0, we assume that there is a cutoff in the wave numbers q 
of the photons; we multiply the density of the states of the electromagnetic fields by e~ q2a2 ^ 4 , 
with a > being a length scale. The cutoff removes all mathematical problems concerning, 
e.g., the exchange of the order of derivatives and integrals, which are abundant in the theory 
without the cutoff. At the end of the calculations we ultimately take the limit a — > 0. 
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In spite of the change in the Hamiltonian and the added cutoff of photon frequencies, 
the analysis of the electromagnetic fields proceeds almost as in ||10|| . In accordance with 
Ref. fI5fl , it emerges from our results that what was called the electric field in [[KJ should 
more properly be interpreted as the electric displacement divided by the permittivity of the 
vacuum e . We henceforth adopt this interpretation. The positive frequency part of the 
electric displacement is expressed in terms of the matter fields as 

D+(rt) = D+(rt) + e f dt' f d 3 r' S(d ge ; r - r', t - t') ^(rV)^(rV) . (4) 

In this approach electric displacement and matter fields are the primary degrees of freedom; 
Dj(rt) is the free electric displacement that would apply if there were no coupling between 
matter and electromagnetic fields. The propagator that takes the radiation from the dipole 
source to the field point is 

S(T>; rt) = — ^- ( d 3 q e -* 2 " 2 / 4 S x ( - x T>) e iqr (e icqt - e~ icqt ) (5a) 
167Heo J q \q J 

° ( P xV)xV ^ r - d| - J ; [r + dl , (5b) 



4vre \l- — i | 

where the delta function has acquired a finite width as a result of the cutoff in the photon 
energy spectrum: 



6 a (x) = exp 



a J 



(6) 



As before, we assume that there is a dominant frequency Q in the problem. In fact, 
we generally assume that a field such as ip e (rt)e' int , and similarly for the electromagnetic 
quantities, varies "slowly" in time in comparison with e~ lQt . From now on a notation such 
as i/) e (rt) refers to the slowly varying field ip e (rt)e int , unless explicitly stated otherwise. 

Based on an implicit cutoff such as a, we argued in |]TD[ that the delta functions with 



plus and minus signs in Eq. ( |5l5| ) conspire to remove a term cx 5(r — r') that results when 
the position derivatives act on |r — r'| _1 . What we did not realize is that this delta function 
does not outright vanish. Instead, it gets smeared to a function whose integral over r is still 



unity but which has a finite width of the order a; see Appendix |A 1| . From now on, as long 



as the integral operator with the kernel S acts on any smooth function cf>(rt)e int of r and 
t in which the exponential is the dominant time dependence, we write 



J d£ J d 3 r' S(£>; r - r', t - t')^{r't')e- iat ' = e~ m J d 3 r' S'(X>; r - r')<f){r't d ) . (7) 

Here 

|r — r'l 

U = t- (8) 



is the retarded time. The monochromatic version of the propagator S, S', may be written 
alternatively 

S'(X>; r) = -L(X> x V) x V — , (9a) 
47re r 

or 

S'(X»;r)=K(X»;r) + Ax» 5(r ). ( 9b ) 

3e 

The final new kernel K(X>; r) is equal to the positive-frequency component of the electric 
field from a monochromatic dipole with the complex amplitude X>, given that the dipole 
resides at origin and the field is observed at r ^ 0. The explicit expression is, of course pi| , 



1 lid™ "1 * / 

K(X>; r) = — {k 2 (hxT>) x n^— + [3n(n • V) - £>](- - ^)e lkr }, (10) 



with 



n = -, k=-. (11) 

r c 

It should be noted that, as the dipole radiation has the 1/r 3 singularity, integrals such as 
(0) are generally not absolutely convergent. According to Appendix [A 1| , we resolve this 
problem by the rule that at least in the immediate vicinity of the divergence r' = r the 
integral is to be done in spherical polar coordinates, and the angles are to be integrated over 
first. It should also be born in mind that the form (]?]) does not apply if the function is 
singular in r. 



We finally consider the quantum expectation value of Eq. (f|). As is always done in this 
paper, we take the free field to be in a coherent state. We also assume that the expectation 
value (Dp) is effectively monochromatic. It is physically evident that, at least in steady 
state, the expectation value of the product V'JV'e will then be monochromatic and a smooth 
function of r as well. Inasmuch as the quantum expectation value of Eq. (|]) is concerned, the 
transformation from the kernel S to the kernel S' is thus allowed. Moreover, the expectation 
value of the free field is a solution to the Helmholtz equation, and the function e lkr /r is 
essentially the Green's function of the Helmholtz differential operator: 

(V 2 + A; 2 )(D+> = 0, (V 2 + k 2 )*—^- = -A-k5(t) . (12) 

In view of Eqs. (|7|) and fl9a]), we thus have from Eq. 

(V 2 + A; 2 )(D + ) = -V x (V x (P+)). (13) 

Classically, the polarization of the medium is defined as dipole moment per atom times the 
density of atoms. It is therefore clear that 

P+(r) = d ge ^(r)^ e (r) (14) 

should be identified as the (positive frequency part of the) quantum mechanical polarization 
operator. 

The value of Eq. (|13"D is twofold. First, it is a local differential equation, as opposed to 
the integral equation (|4]). Second, it has a well-known counterpart in the classical electro- 
dynamics of polarizable media. This reinforces the interpretations of D and P as electric 
displacement and polarization operators. 

3. Matter fields 

We now turn to the equations of motion of the matter fields. Under the assumptions that 
the density of excited atoms is low and that an atom moves much less than a wavelength of 



S 



resonant light during the time it remains excited, we have the equations of motion for the 
fields describing excited and ground state atoms, 

ij> e (rt) = i6 ^e(rt) + | d eg ■ E+(rt) ^(rf) , (15a) 



fart) = iE"(rt) -&Mrt) -<S^£) ^( rt ) + | 



^(rt) • (15b) 

c 



As the notation implies, 



E + (rt) = - [D+(rt) - P + (rt)] (16) 

is to be interpreted as the electric field. Furthermore, 5 = Q — loq is the detuning of the 
characteristic frequency of the light Q from the atomic resonance frequency Uq. We have 
carried out the Rotating Wave Approximation that takes into account the dominant field 
frequency Q. Finally, H c m is the one-particle Hamiltonian governing the cm. motion of the 
atoms, and the time derivative with the subscript C represents collisions. 



4- Summary remark 

The two-state model of quantum optics is immediately seen to underlie Eqs. (|T5|), and 
Eq. ([|) describes the total field as the incident field plus the fields radiated by the dipoles, 
complete with propagation delays. In spite of the familiar appearances, though, it should 
be noted that the only real approximations so far have been to ignore the cm. motion and 
collisions of the excited atoms. The formulation still fully accounts for quantum statistics 
of the many-atom system, and for the quantized electromagnetic fields. The effect of the 
dipole-dipole interactions on the transition frequencies and linewidths of the atoms is still 
included. On the other hand, as we have ignored the cm. Hamiltonian for excited atoms, 
collisions between ground state atoms and excited atoms can no longer be discussed. 
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B. Eliminating the vacuum field 



Even in the absence of applied electromagnetic fields the atoms bathe in vacuum fields 
that cause spontaneous emission and Lamb shifts. The purpose of the present section is to 
account for the vacuum fields. While pursuing this goal, we need to be prepared for singular 
functions with rapid spatial and temporal variations. We therefore start with the general 
field equation (H). Moreover, for the time being we argue in terms of the original atomic 
and electromagnetic fields untempered by the exponential e tnt . 

To begin with, we insert Eq. (ffl) into Eq. (|15a|), and obtain 



ip e (rt) = -iu ip e (rt) 



+ ^-d eg • | D+(rt)^(rt) - d 9 v^(rt)^(rt)Vv(rt) 
+ e £ dt J d 3 r' S(d gV ; r - r', t - t')^ gl {v't!)^ el {v't')^ 9 {rt) } . (17) 

We are ultimately interested in quantum expectation values of atomic and electromag- 
netic field operators, and thus wish to be able to take expectation values of expressions such 
as Eq. (|lT| ) easily. It would be especially valuable to have the free- field operator D^(rt) 
farthest to the right. Because the (initial) free electromagnetic field is assumed to be in the 
coherent state, in an expectation value this operator would then reduce to a multiplicative 
classical field: the relation 

(ODj(rt)> = <0)<D+(rt)> (18) 

applies to any operator O. Evidently we need commutators between atom operators and 
free-field operators, so that we may move the latter to the desired positions. 

Any atom operator, of course, commutes with the total electric displacement operator 
D ± at the same time. We thus have from Eq. (|]) 

Q=[V s (rt),D+(rt)] 

= -e f dt 1 /rfVS(d 3 , e ;r-r',t-t / )[^(ft),^(r , t , )^ e (r / t / )]. (19) 

J —CO J 

10 



Here we are preparing for the eventuality that the commutator is required for two different 
field points. 

The standard way of dealing with vacuum fields in quantum optics is the duo of Born 
and Markov approximations: the atom operators evolve as if no electromagnetic fields were 
present (Born) during the short vacuum correlation time (Markov) |15|| . For implementations 
of this idea in the Heisenberg picture see, e.g., [|I^ -[l8||. We evaluate the commutator Q under 
an approximation which, we think, is the field theory equivalent of the standard Born and 
Markov approximations: We assume that during the time it takes radiation reaction effects 
to assert themselves, the atom fields evolve as if they were completely noninteracting. We 
temporarily restore the explicit notation for magnetic quantum numbers, and write the 
Born-Markov approximation for Eq. ([19|) as 

Vw(rO = 4f? E e^*-*^ b gmk (t), VeM(rt') = £ e^-^'-^b^t) . (20) 
vv k W k 

Here ek = hk 2 /2m gives the dispersion relation for an atom with mass m, b are boson 
operators, and the sums run over the wave vectors k appropriate for the quantization volume 
V. In the standard continuum limit the relevant commutator becomes 

& 9m (HUl m ,,(rn = J^l d " k c*-^-*^) . (21) 

We use S from (|5a|), the commutator from (^1|), ijr gm i from (p0[), and add the conventional 
convergence factor e~ vt to the time integral. Equation (O) is cast into the form 

Q = ]T e &eM'K / dre-v / d 3 q 

167rV^ jttd' Jo J 



e iq-(r-r) g-aV/4 Q x j S x d mM , ) (e^ T - e"^) e^-^+^r _ fry 

q \q J 

Our final approximation is to ignore the cm. energies in comparison with the energy of 
the internal transition of the atom; we write ujq — ck q — — ^- The results is interesting: 

roc 

Q = -eJ dr S(d 9e ; f - r, r)e inT ^(rt) . (23) 
Jo 

The time integral is the same as the definition of the kernel S'(d 9e ; f — r) in (|9[), albeit still 
containing the cutoff parameter a. The cutoff is truly needed: in our immediate application 
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to Eq. (|17]) we are to set r = r, and without the cutoff we would have to contend with a 
pernicious singularity of the type S(r)/r 3 . For a small but nonzero a, the result is 

[Vfc(rt), Dj(rt)] = d ge - ^(rt) • (24) 

Given the sum rule for the dipole moment matrix elements, ( |A15|) , the relevant terms in 
Eq. ( |17|) become 



4> e (rt) = -iu ip e (rt) + — d eg ■ Dj(rt) ^ fl (r£) 

/2Co 



-iw ^ e (rt) 



7 + *; 



i/j e + ^-d eg -i/j g (rt)B+(Tt)... , (25) 



37T 2 eo^.a 3 

where X> is the reduced dipole moment matrix element. The imaginary part in the second 
term on the right-hand side of Eq. fl25|) diverges as the photon momentum cutoff goes to 
infinity with a — > 0. This part, after a proper renormalization, contributes to the Lamb shift. 
From now on we assume that the Lamb shift is already included in the transition frequency, 
and ignore the a~ 3 term in fl2"5p. What remains is the familiar spontaneous linewidth of the 
atomic transition, 



V'co 3 

,3 



7 = ^T^- ^ 



There are no divergence problems with the commutator of ip g and if the position 
arguments are different. We simply write 

[^(rf), D+(rt)] = -e S'(d 9e ; r - f)^ e (rt) . (27a) 

In fact, if the divergent in-phase part of the dipole field at f = r is ignored (or incorporated 
in the Lamb shift), we may interpret Eq. (|27aj ) to be valid even for f = r. 

By the same token, we may carry out all commutators between atom fields and free 
electromagnetic fields. The two additional nonvanishing commutators that play some role 
in this paper are 

[^t(ft),Di(rt)] = eo ^J(rt)[S / (d ae ;r-f)]*, (27b) 
[^ e (rt),D F (ri)] = -e i/; g (H)[S'(d ge ;r-r)}*. (27c) 
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At this point we restore our convention of slowly varying fields. We also add another 
assumption to the effect that light has ample time to propagate across the atomic sample 
during the time that the slowly varying fields change appreciably. This permits us to ignore 
propagation delays in the time arguments of the slowly varying fields. The retarded time td 
is simply replaced by the external time t. We have thus obtained an equation of motion for 
the excited state atom field that contains an explicit radiative damping, 

ij) e (r) = (iS - 7)V> e (r) 

+ ITo ^ ' { ^ r ) D ^ ~ d 9' e '^'( r )^ e '( r )^( r ) 

+ 6 / rfV S'(d g v; r - r')4(r> e ,(r')^(r) j . ( 28 ) 

Here, and in our subsequent expressions, the common time t is omitted in the notation. 

Unlike in the ordinary treatments of spontaneous emission in the quantum optics of an 
isolated atom, no short vacuum correlation time suggests itself in our formulation. The 
use of the free-field evolution as in (|20|) during the "vacuum correlation time" may thus 
seem like an ad-hoc assumption. This approximation, however, did produce spontaneous 
damping and Lamb shift in accordance with the one-atom theory. Of course, even in standard 
quantum optics the atomic variables do not evolve completely freely during the vacuum 
correlation time. Spontaneous emission itself, as well as external driving electromagnetic 
fields and collisions between the atoms in principle affect spontaneous emission, but at 
ordinary conditions for laser spectroscopy these influences are negligible We conjecture 



that the same applies in our field theory version of spontaneous emission. Finally, ignoring 
cm. energies in comparison with the energy of the atomic transition is nothing new either. 
This is a standard approximation in the derivation of spontaneous emission in the theory 
of light pressure If such energies were included, a velocity dependent spontaneous 

emission rate would emerge in manifest contradiction with special relativity ||20|| . 
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C. A hierarchy for operator products 



A particularly relevant atomic variable is the polarization operator, (14), which acts as 



the source for secondary radiation. Generalizing slightly, we now embark on a study of the 
time evolution of the operator product ^(r)?/' e (r). 

We have in mind situations in which collisions and cm. motion of the ground state atoms 
have come to a steady state before the driving light is turned on. We regard the external 
field as a small perturbation, so that the ground state atoms remain materially unperturbed 
in the presence of the light. As the final item, we assume that collisions and cm. motion 
of the ground state atoms take place on a time scale much longer than the spontaneous 
emission time scale 7 _1 . They are therefore not expected to interfere with spectroscopic 
probing of the atomic transition. Under these assumptions we will henceforth ignore the 
collision terms and the cm. evolution of the ground state atoms altogether. Nonetheless, 
mathematical consistency dictates that in the intermediate steps of our calculations we take 
into account some light driven evolution even for the ground state atoms. 

We thus have the equation of motion from Eqs. ( |15b| ) and (f28|), 



d eg , • { VjVyD+ - ^JP + ^' + e / d 3 r>t S '(P+( r '), r - r')Vy 
~i de ' 9 ' { V^ e D£ - ^P + ^e + eo J dV4S'(P + (r'),r-r')i } , (29) 

where we have shown explicitly only the nonlocal position dependence. Generalizing, we are 
evidently about to derive a hierarchy of equations for operators of the type P + , ^P + -^ 9 , 
^t/^P 4 'ip gip g , . . . , with different position arguments for the different fields. 

Now, light has to be present in order to produce excited atoms. Each excited state 
field corresponds to one order in the strength of the driving light. To first order in we 
might thus ignore the second term on the right-hand side of Eq. ( PP[ ) altogether. This is, 
generally speaking, what we will do: we only retain those products of operators that involve 
at most one of the operators Dp, D^, ip e or However, caution must be exercised for two 
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reasons. First, we will eventually arrange all atom fields to normal order; creation operators 
to the left, annihilation operators to the right. Besides, we move the free-field operators 
to prescribed positions as well. In the process commutators are generated that may be of 
different order in the strength of the driving field than the original terms. Second, some of 
the commutators are flat out divergent, analogous to the Lamb shift. It may be shown that 
such extra Lamb shifts cancel exactly, order by order in the strength of the driving field, but 
the cancellation of course fails if the calculations are not consistent in the orders. We will 
not dwell on the latter aspect anymore, but simply ignore all orders higher than the first 
one immediately at the point when the operators have been brought to the desired order. 

We illustrate the process of deriving the hierarchy of equations for operator products 
with a detailed treatment of the time evolution of a particular product; 

d 



dt 



Vi(r>J(r)^(r)^(r') 
d 



dt 



V£(r)V' e (r)Vy i (r')+V' 



4(0 



di 



^(r)Ve(r) 



1MO 



d 



+ i(r>J(r)i(r)^;(r'). 



(30) 



In the way of preparation, let us first note from Eq. @ and its hermitian conjugate that 
the electric displacement may be written in the form 



D ± = D± + Df , 



(31) 



where the source term is a normal-ordered combination of atom fields. We begin our 



analysis with the third term on the right-hand side of Eq. (p0|). By virtue of Eqs. ( |15b|) 
and ( |I6| ) we first have 



d 



dt 



^(rO^( r )^e(r)^.E-(rO^(rO 



— ^(rO^(r)ik(r)d 



9[e' 



D-(r') - p-fr' 



he 



^ gi (r')^ g (r)d g[e , • D-(r')^ e (r)^(r / ) - ^(^(r^tod^, • P-(r')Vv(r' 
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<(r')^(r)dv e ,-D^(r')^ e (r)^(r') 



i 

+ Vi(r')V>J(r)d gie , ■ D;(r')^(r)Ve'(r') - V£ Wjto^tt'W ' P-(r')Vv(r')] • (32) 

Since the total displacement D~ commutes with all atom fields, we have first moved it 
between atom creation and annihilation operators. The term involving the source field 
is then readily in normal order. Besides, here the source field term is third order in the 
perturbation, so it may be omitted. However, we are not yet done with operator orderings. 
First, just as the free-field operator is profitably moved to the right of atom operators, the 
free-field operator should be transported all the way to the left. By virtue of Eq. (|27b[) . 



this leaves behind two commutator terms. Nevertheless, both the term with remaining 
and the commutators are formally third order in the strength of the driving field, and we 
ignore them all. Finally, the term involving P~ is not yet in normal order, but it may be 
made so easily by using the commutators of the atom fields. The rearranged term is third 
order and negligible, but the generated commutator is first order. We eventually have 

Vi(r>J(r)Var)^(r') 



~ d g , e , ■ d eg „5(r - r')Vj 1 (r')^(r)V 9 "(r / )V'e'(r / ) . (33) 



dt 

i 

he 

The same analysis may be carried out with the other two terms in Eq. fl3~0|). The first 
term contributes nothing in the first order, while the second term gives a homogeneous term 
proportional to iS — 7, a driven term proportional to Dj, and something of a two-atom 
analog of radiation reaction. The final result is 
d 



dt 



^ i(r >t (r) ^ (r) ^, (r0 

[iS - 7 )^ 1 (r')^( r )^(r)^(r') + ±- {^(r>J(r)Vv(r)^( r, )<W ' Dj(r) 



+4(r')^(r)V 9 '(r)^(r')d efl ' • e S'(d g ur, r - r') - d g[e ,8(r - r') 



hc 

-^ gi (r')^ g (r)d egl • P+(r)Vy(r)^(r') 

+e ^ 1 (r')^( r ) /^V'd e ^-S'(P + (r");r"-r)^(r)^(r')} . (34) 
Continuing in this manner, we obtain the equations of motion for an entire hierarchy of 
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products of atomic operators. As before, we put the positive-frequency free-field operator 
to the right, all atom operators to normal order, and then keep only the terms that are first 
order in the perturbation strength. The full result is 

= (iS - 7)^1 (r„) • • • ^ 1 (r 1 )V e (r 1 )^/(r 2 ) . . .^(r„) 



+J— {^lS r n) ■ ■■^l 1 (r 1 )^ g (r 1 )^ 2 (r 2 ) . . .^(r n )d e5 • D£(ri) 

+e Vi(iv) • • . ^(ri) / d s r'd eg ■ S'(P+(rVi - r')^( ri )^(r 2 ) . . (r n ) 

n 

+e ipl n (r n ) . . .Vi(ri)V a (ra) ^ ^(r 2 ) . . .^^{T^^Y^^r^) . . .^;(r n ) 

fc=2 

xd ej -W(d 9 ^;ri-r ]fe )} , (35) 

where the notation in the last term implies that ip g > (r^) is missing from the /c term on the 
sum. We have defined 

W(X>; r) = S'(X>; r) - - T> 5{v) 

= K{V;v)-^-'D8{r). (36) 

This is precisely the classical expression of the electric field (not displacement) of a dipole 
T> residing at the origin, as measured at the point r. Even the peculiar delta function 
divergence of the dipolar field at the origin fbl is there. 



D. A hierarchy for correlation functions 



By taking expectation values of the operator hierarchy fl3"o|), we obtain a hierarchy of 
equations of motion for correlation functions. In order to simplify, in the rest of the paper 
we only consider a J g = — > J e = 1 transition. Then there are no Zeeman substates in the 
ground level, and a single g suffices in all of the Eqs. (^). The three excited Zeeman states 



are also handled easily, c.f. Appendix |A 2 . 
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We define a succession of correlation functions 



P 3 (r 1 ,r2;r 3 ) 



P 2 (ri;r 2 ) 



(V»2(Pi)d^e(ri)> = (P + (ri)) , 
(^(r 1 )P+(r 2 )^(r 1 )} , 
(^(r 1 )^(r 2 )P+(r 3 )V s (r 2 )V fl (r 1 )} , 



(37) 



and similarly 



p 2 (ri,r 2 ) = (^(11)^(12)^^2)^^1)), 



(38) 



Pjfc(ri, . . . , rfe_i; rfc) is the correlation function of polarization at and atom density at 
k — 1 positions r 1; . . . , r^-i, and pk is a fc-point density correlation function. All of these are 
normally ordered. 

We finally reiterate that the driving field is in a coherent state, so that the factorization 
(|i~8|) is warranted. In fact, without further ado, we let stand for the expectation value 
of the coherent free field, or, equally well, for a classical incident field. It is now a simple 
matter to derive a hierarchy of equations for the correlation functions from the operator 
equations (j35|) . We consolidate the results into the form 



P 3 (ri, r 2 ; r 3 ) = (iS - 7)Ps(ri, r 2 ; r 3 ) + G(r 3 - ri)P 3 (r 3 , r 2 ; n) + G(r 3 - r 2 )P 3 (ri, r 3 ; r 2 ) 
+mp 3 (ri,r 2 ,r 3 )D£(r 3 ) 



Pi(; n) = (i6 - t)Pi(; t x ) + i/«p 1 (r 1 )D+(r 1 ) + J Sr 2 G(ri - r 2 )P 2 (n; r 2 ) , (39a) 
P2(ri; r 2 ) = (id - 7)P 2 (r : ; r 2 ) + G(r 2 - ri)P 2 (r 2 ; r x ) 



+mp 2 (ri,r 2 )D^(r 2 ) 




(39b) 




(39c) 
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We have defined the scalar constant 



V 2 

k = —, (40a) 
he 



and the 3x3 tensor 



Gij(r) = in 



99 S ^2 



gifcr 



47rr 



M( r ) . (40b) 



9rj 9r 3 - 

The first two equation of the hierarchy ([39]) coincide exactly with the results of Morice 
et al. ||. These authors do not proceed any further, but their method, which at this point 
had become tantamount to classical electrodynamics, undoubtedly could have yielded the 
entire hierarchy. 



The terms in (|3~9 ) with id — 7 obviously come from the damped free evolution of the 
polarization in each correlation function, and the term oc Dj, corresponds to excitation of a 
ground state atom by the driving light to make polarization. To grasp the dual role of the 
tensor G, let us consider the equation of motion for P 2 (r!;r 2 ) as an example; correlation 
function of polarization at r 2 and density at ri. The integral term obviously characterizes 
processes in which a dipole at yet another position r 3 radiates and thereby promotes a 
ground state atom at r 2 , so that density becomes dipole density. On the other hand, the 
term with G(r 2 — ri) describes photon exchange between the two sites ri and r 2 ; an excited 
atom radiates at ri and falls to the ground state, while the emitted radiation promotes an 
atom at r 2 to the excited state. 

We have implemented several approximations. The most relevant physical assumption is 
the perturbative limit with respect to the strength of the driving light, the most conspicuous 
technical assumption is the J g = — > J e = 1 transition. However, Eqs. (|39|) do not contain 
any assumptions concerning spontaneous emission except for our field theory version of 
the Born and Markov approximations. Moreover, there are no assumptions, ad hoc or 
otherwise, concerning multiple scattering of light or resonant dipole-dipole interactions; these 
are included exactly. Hence, so are collective linewidths and line shifts. 

As far as the interactions of atoms with electromagnetic fields are concerned, we have 
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regarded the atoms as point dipoles. For real atoms at short distances, when higher multi- 
poles and electron exchange become relevant, this assumption evidently fails. The relevance 
of contact interactions and 5 function contributions to the dipolar field is questionable, as 
both operate at zero distance between the atoms only. Now, real atoms cannot overlap 
because of the hard core of the interatomic potential. A reader troubled by the delta func- 
tions may therefore want to consider cutting off and setting to zero all correlation functions 
at distances between the atoms shorter than the typical length scale of a molecular bond. 
The effect is that all delta function contributions to the field propagator G of (|40b| ) should 
be omitted. As the derivatives of e lkr /r also produce delta functions, such an omission is 
tantamount to replacing the term —5ij5(r) by — |5jj5(r). However, in the present paper we 
use the propagator exactly as given in Eq. (|40b|) . 

III. EXAMPLES 

In this section we illustrate the correlation function hierarchy with a few simple examples. 
At this time we have made little progress toward a full, exact solution of the hierarchy ( |39| ) 
in any nontrivial situation. Evidently, radically new approaches are needed. We hope that 
either ourselves or our readers will in the end be inspired to come up with a successful 
solution of, say, the optical response of a Bose-Einstein condensate (BEC) in the limit of 
truly dense sample, pA 3 1. 

A. Semi-infinite BEC without collective coupling 

1. Nature of condensate 

An ideal Bose condensate of noninteracting particles is made of a macroscopic number of 
particles in the same one-particle quantum state. Traditionally, the condensate is described 
by a macroscopic wave function, whose absolute square gives the particle density and which 
also has a phase. More in the vein of quantum optics, one could assume that the condensate is 
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in a coherent state, albeit with an unknown phase. In normally ordered operator expressions 
the field operator tp g then behaves as a c-number. We write 

p Jfc (r 1 ,...,r Jfc ) = <^(ri)...V,(ri)> (41) 

~ ^*(n) . . . V*(r fc )V s (r fc ) . . . i) g {vi) = V*(ri)V s (ri) . . . ip* g (r k )ip g (r k ) 

= (Vj(ri)^(ri)) • • • (^(r fe )^(r fe )) = Pl ( ri ) . . . Pl (r fc ) . (42) 

In other words, density correlation functions factorize. For the purposes of the present paper, 
we take the factorization of normally ordered density correlation functions as the hallmark 
of the condensate even for a weakly interacting condensate. We ignore noncondensate atoms 
altogether. 

Contrary to the experimental realities, we ignore the finite dimensions of the condensate. 
We cannot outright declare pi(r) as a constant all over space, because propagation through 
an infinite medium would cause extinction of all light before it reaches any position with 
finite |r|. Instead, we assume that a homogeneous condensate with density p fills the half- 
space z > 0. We assume that the incident and induced radiations as well as the induced 
polarizations all propagate in the z direction. Finally, as we have the m g = spherically 
symmetric ground state that cannot exhibit any directional preferences, we take all fields to 
have the same transverse polarization e. In particular, the initial free field is written 

D+(r) =D F ee ikz . (43) 

This has the dispersion relation of light in vacuum, an oddity in the presence of matter. In 
fact, in accordance with the Ewald-Oseen extinction theorem |2l|], it will turn out that the 



matter responds with a field that exactly cancels the applied field. 

2. Optical response without collective coupling 

We solve the response by ignoring the collective line shifts and dampings, i.e., those 
G terms in Eqs. Q59| ) that do not appear inside integrals. It turns out that the simplest 
conceivable ansatz, a fully factorized, damped plane wave solution of the form 
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Ppr-iee**", z 1 >0,...z n >0; 
P n (ri, . . . ; r n ) = { (44) 

0, otherwise. 

succeeds for suitable choices of P and k', with $$(k') > 0. 
To see this, we need the integral 



/ Sr 2 e* ■Q{r 1 -v 2 )-ee ik ' Z2 =m 

Jz?>0 



k' 2 



k' 2 -k 2 2k(k-k') 



(45) 



where the vector e takes care of the polarizations. The integral is valid for Z\ > 0. For 
z\ < the integral yields a reflected wave oc e~ lkz instead, but we do not consider this case 
any further. With the ansatz ([14]) , the integrals on the right hand sides of (|39p produce 
sums of two exponentials, one with the wave number k appropriate for light in vacuum and 
the other with the wave number k' for light in the medium. In steady state of Eqs. fl39|) the 
vacuum component oc e lkz must cancel the corresponding free-field terms, and the remaining 
e ik'z t erm mus t pair up with the polarization correlation functions P n . All of the Eqs. ( |5UD 
then reduce to the following two conditions: 

{ iS -i + &) p=0 ' (46a) 

D ^2W^) P = °- (46b) 
The first one gives the wave number k', and the second one may then be read as a condition 
for the polarization amplitude P. 

One expects that the total displacement from (|J) should be of the form D e e tk ' z as well. 
With the choice (|46b|) , the vacuum type contributions e lkz indeed cancel. We have the 
condition for the polarization amplitude and the amplitude of electric displacement, 

D -w^ p - < 47 » 

This implies that an electric field of the form E e e lk ' z also propagates in the medium, with 
the amplitude given by 

B =^f = I^P. ,48, 
e e k — k z 
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As customary, we define refractive index n in such a way that k' = nk, and susceptibility x 
such that P = e xE. Equations (EH)-(E51) immediately give the explicit expressions 



» 2 -1 = X = -t^t- (49) 

for these quantities. 

The remarkable feature of the result ([49f ) is that it is so unremarkable: susceptibility 
is obtained as atomic polarizability times atom density. This is precisely the conventional 
column density approach that the experimenters routinely use to analyze their BEC results. 
In other words, we have proven that, for a plausible model of the BEC and within a precisely 
formulated approximation that ignores collective linewidths and shifts, the column density 
approach is exact. We regard this as an important result, in that it displays precisely and 
explicitly the underlying assumptions of the column density arguments. The flip side is that 
at high density simply ignoring the collective effects is another uncontrolled approximation. 

B. Density expansion to second order 

In the present section we review the density expansion of Morice et al. || from the point 
of view of our development. In effect, they truncate the hierarchy for correlation functions 
by writing 

P 8 (n, r 2 ; r 3 ) ~ P -^A rjJ ) . (50 ) 

If (and probably in some fairly strong sense only if) one inserts this particular factorization 
into (|39b|) , one may eliminate at the same time both the free- field term and the integral 
term from the steady-state versions of Eqs. (|39aj) and ( |39 b|) . This gives an algebraic relation 
between Pi and P2, 

= P 2 (r i; r 2 ) - 5fep*) P 2 ( r2; Pl) + E^l P l(; r 2 ) . (51) 
{16 -j) pi(r 2 ) 

Now assume a constant density of atoms p in the half-space z > 0, and write the 
physically justifiable ansatz 
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Mri, r 2 ) = p 2 [1 + <p(ri - r 2 )] . 



(52) 



The relation (^Tj) and its counterpart with 1*1 and r 2 interchanged may then be solved to 



give 



P 2 (ri;r 2 ) 



p[l + fjjcx - r 2 )] 
1 - [G(n - r a )/(i* - 7)] 1 



PiO^-^^PxOn) 

(zd - 7) 



pPi(;r 2 ) 



[(^ 2 )P 1 (;r 2 )- 5 P 1 (;r 1 )] 



(53) 



1-^2 

where we have introduced an obvious temporary notation. The reason for the split in the 
second form of ( |53| ) is the following. Suppose we attempt an ansatz of the form (fi3D, 
together with Pi(; r) = P ee lk ' z . Then, when (|53|) is inserted into the integral in right-hand 
side of (|39a|) , in accordance with ([45 ) the term pPi will produce a term with the spatial 
dependence of the free field e lkz . All other contributions to the integral will behave as e tk ' z , 
at least for z — > oo. Now, the requirement that the free-field contributions cancel (Ewald- 



Oseen extinction theorem [21] again) gives one relation between Dp, P and k', and Eq. (p~3|) 
furnishes another. One then obtains an equation out of which one may solve k', and therefore 
ultimately the entire response of (a thick slab) of the gas. This equation reads 

1 



k 2 



5 + 1 + C 



(54a) 



with 



C 



P 



d 3 r V (r)e- tk ' z e* • G(r) ■ e 

- ik ' z [Q{v)/{i5 - 7 )] 3 - [G(r)/(i5 - 7 )] 2 



iS — 7 

+p [ d 3 r [1 + <p(r)] e* 



e . 



(54b) 



l-[G(r)/(<tf- 7 )] 2 

Ingenuous as the analysis of Ref. P is, it elicits two questions. First, the factorization 
(p0|) is the way it is to facilitate the mathematics, not because ( p0|) would be a particularly 
apt physics assumption. There is no guarantee that the correlations between the ground 
state atoms are treated adequately. Second, while Morice et al. argue that the result is 
the correct expansion in the parameter pA 3 up to the order (pA 3 ) 2 and, in fact, includes 
all multiple scattering events between any pair of atoms, the mathematical structure of the 
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hierarchy fl39|) does not directly bear this out. The difficulty is that each integral of the form 
/ rf 3 r/ c GPfc must produce a component that cancels the associated free-field term. In other 
words, the integral does not automatically signal an increasing power in pA 3 . Within the 
present approach, a rigorous mathematical counting of the powers of pA 3 seems elusive. 

C. Resonant dipole-dipole interaction of two atoms 

Let us take two atoms with small non-overlapping cm. wave packets 0±(r) centered 
around r±. Inasmuch as the sizes of the wave packets are much smaller than the wavelength 
of the exciting light, for the purposes of the analysis of optical response we may write the 
atom density and the density correlation function as 

Pl (n) = |0 + ( ri )| 2 + |0_( ri )| 2 ~ 5( ri - r+) + 6{n - r_) , 

P2 (r 1 ,r 2 ) = |0 + (r 1 )ri0-(r 2 )r + |0 + (r 2 )| 2 |0-(ri)r 

~ 5(ri - r + )5(r 2 - r_) + 5(r 2 - r + )5(ri - r_) . (55) 

To simplify the results further we assume that the geometry of the situation is such that 
for the driving light at the positions of the atoms we have D^(r + ) = D^(r_) = Dp. In 
particular, this is justified if the atoms are well within a wavelength of one another. In such 
a case the two-atom system may be discussed in terms of Dicke states, some superradiant 
and some subradiant. The field configuration we have chosen does not excite the subradiant 
states at all, so these will not come up in our analysis. 

Since we have two atoms present, all correlation functions referring to more than two 

atoms vanish: p n = and P n = for n = 3, 4, The steady state of the hierarchy (|39| ) is 

then found trivially: 

PxO rj) = p[6( ri - r+) + <f(n - r_)] , (56) 

with 
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Unraveling Eq. fl57|) (with the 3x3 tensor G in the denominator) gives a complicated 
expression that depends on the polarization of the driving light Dj, direction between the 
atoms, and the distance between the atoms r. The main features, however, are intuitive 
obvious. For r C A the optical response shows a split resonance at 

5 =-7wT- (58) 



47rr 3 27rr s 

The configurations that produce only one or the other of these resonance are such that 
the polarization of the driving field (and hence, of the induced dipoles) is perpendicular 
or parallel to the vector joining the dipoles, respectively. These correspond the doubly 
degenerate tt u and the nondegenerate a u configurations of the molecule consisting of the two 
atoms. The cooperative width of the resonances is 27 instead of the one-atom linewidth 7. 
Incidentally, the subradiant states that are invisible within our approximations correspond 
to the molecular configurations n g and a g , which appropriately do not have dipole coupling 
to the molecular ground state. 

Because our theory is linear in the external field, it should not come as a surprise that 
exactly the same dipole-dipole response ensues for two classical, isotropic, charged harmonic 
oscillators. In such a calculation one may want to put in the one-oscillator damping in 7 
by hand in order to dispense with an explicit treatment of the divergent radiation reaction. 
Nonetheless, the collective linewidth and dipole-dipole interactions are easily derived from 
the classical analysis of the radiation that the oscillators exchange among themselves. On the 
other hand, the molecular analogy also hints at some aspects of physics that are missing from 
our formulation. First, once more, at close enough distances two atoms may no longer be 
regarded as point dipoles. Second, since we have ignored all cm. evolution, the possibility of 
quantized vibrational states of the molecules and the ensuing optical resonances have fallen 
by the wayside. 

As the connection of our hierarchy to dipole-dipole interactions has now come up explic- 
itly, we should point out that there are numerous discussion of two- and n-atom response 
in the literature in which (in effect) the dipole-dipole interaction is derived (in effect) by 
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eliminating vacuum electromagnetic fields with the aid of Born and Markov approximations; 



see, e.g., |16[. As it comes to dense and/or degenerate gases, it seems that dipole-dipole 
interactions are going to occupy an increasingly prominent position in theory and eventually 



perhaps in experiments as well. "Nonlinear atom optics" [p2} -p5[|, in which the self-interaction 
of the atom wave derives from dipole-dipole interactions, is a prominent example. In fact, 
the analysis of Zhang and Walls [^3j in terms of Heisenberg picture field operators is in spirit 
quite close to our treatment, and it seems plausible that our correlation function hierarchy 
could also be derived from the Schrodinger picture master equation as presented by Lenz et 
al El. 



D. Prospects for exact solution 

Dipole radiation presents mathematical difficulties both at short and long distances. One 
has the 1/r 3 and indeed a delta function divergence at short distances, which means that 
the results are sensitive to short range correlations between the atoms. Also, the dipole 
interaction falls off as 1/r at large distances. Integrals involving the dipole interaction are 
not absolutely convergent on the fall-off of the interaction alone, and local approximations 
of the type 

J dVG(r - r')/(r') ~ /(r) J dVG(r - r') 

cannot be made. There may be a global coupling between the electromagnetic fields and 
polarization reaching across the entire sample. At the face of such mathematical hazards, 
any uncontrolled approximation should be viewed with suspicion. An essentially exact, most 
likely numerical, solution of the hierarchy fl39|) appears highly desirable. 

The idea of attempting a direct solution of the hierarchy of integral equations (|3^) in any 
conceivable future is clearly stillborn. Instead, on several occasions we have seen hints for 
an alternative. We have noted in our treatment of two atoms that in the limit of low light 
intensity the atoms act like classical charged harmonic oscillators. Moreover, even though 



27 



the derivation was in terms of commutator properties of various fields, both the shift and the 
damping of the atoms associated with the dipole-dipole interaction could be viewed simply 
as manifestations of the classical radiation transmitted from one atom to the other. It is 
conceivable that a proper (stochastic) spatial distribution of classical radiators could share 
the correlation function hierarchy fl5Pp. The hierarchy might thus be solved by simulating 
a system of classical atoms and classical electromagnetic fields numerically. However, so far 
we have no mathematically rigorous prescription for such a simulation, let alone a practical 
implementation. 

IV. CONCLUDING REMARKS 

We have presented a fully quantum mechanical, careful analysis of the response of a gas, 
possibly degenerate, to electromagnetic fields. The main technical ingredients are the field 
theory version of Born and Markov approximations, and procedures to move atom fields and 
electromagnetic fields to a certain (basically normal) order. The outcome is a hierarchy of 
equations of motion for atomic correlation functions, in this paper specifically developed for 
the limit of low light intensity. 

Under our assumptions, notably low intensity, the hierarchy is in a sense obvious. In 
retrospect, it could have been outright guessed on the basis of classical electrodynamics. 
We note, though, that our methods would work in many generalizations that go beyond 
classical physics. Given that even the simplest low-intensity limit has not yet been solved 
satisfactorily for near resonance response of a dense gas, we do not address more complicated 
cases at any length. Nonetheless, a few possible generalizations should be mentioned. 

The price for an arbitrary intensity would be a more complicated hierarchy containing 
atomic correlation functions with more than one excited state field ip e . In the case of an 
arbitrary intensity, a two-level atom no longer behaves identically to a charged harmonic os- 
cillator, scattered light may be nonclassical, and statistics of the atoms may play a nontrivial 
role. The hierarchy for an arbitrary intensity most likely will not admit a classical simula- 
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tion. Solutions for a finite number of atoms and density expansions might be extracted, but 
the full hierarchy would present a truly daunting problem. 

While we have treated the cm. degrees of freedom quantum mechanically and at all times 
properly retained the quantum statistics of the atoms, we have for the most part ignored 
the cm. Hamiltonian. We have effectively consigned the atoms to immobility. This is by 
no means necessary. We could add kinetic energy of the atoms, a confining potential, and 
even molecular potential curves to the theory. Of course, this again entails complications: 
the entire physics of molecules made of a single atomic species, photon recoil, cold collisions, 
etc., become special cases of our approach. 

The ultimate objective of our formulation is to find the expectation value of the polar- 
ization, (P + ). Out of (P + ) one may deduce the expectation value of the scattered field, 
and hence, the expectation value of the total electric field (E + ). The flaw here is that a 
typical detector of light does not measure (E + ), but rather expectation values of quadratic 
quantities such as (E~E + ). One may approximate, say, the measured intensity as 

I(rt) = (E"(rt)E + (rt)) ~ (E~ (rt)) (E+ (rt)) , (59) 

but it is known in quantum optics that this type of an approximation broadly speaking misses 
the intensity of inelastically scattered light. Our hierarchy is tantamount to a collection of 
classical linear harmonic oscillators interacting with light, a system in which one expects 
elastic scattering only. However, when one goes beyond the low-intensity limit, even a single 
two-level atom scatters inelastically; and if the cm. hamiltonian is fully included, photon 
recoil gives additional inelastic scattering. To include inelastic scattering properly, one needs 
to develop a hierarchy starting from (P"(rt)P + (r't)), and proceed consistently at least in 
the second order in the strength of the driving electric field. 

Finally, to calculate the spectrum of scattered radiation, one employs the two-time cor- 
relation function of the electric field (E~(r£)E + (ri')), which is obtainable from the two-time 
polarization correlation function (P~(rt)P + (r't / )). To compute the latter, one needs not 
only the field theory version of Born and Markov approximations, but also further consid- 
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erations that essentially amount to the regression theorem [EH . 



The modifications of the spectrum of the scattered light in the presence of the condensate 



discussed in Ref. |10| are due to inelastic scattering associated with photon recoil. Ironically, 
to analyze even this seemingly simple case within our present framework, we would have 
to generalize so that the cm. motion is included, and also develop the quantum regression 
theorem. Our linear hierarchy is no panacea; major generalization are needed in many 
relevant problems. Nonetheless, we hope that the eventual solutions of the hierarchy, by 
means of classical simulations or otherwise, will shed new light on the near-resonance optical 
response of dense atomic samples. 
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APPENDIX A: MATHEMATICAL DETAILS 

1. Dipole radiation with wave number cutoff 

We first consider the following integral involving the propagator that governs dipole 
radiation from matter fields: 

/•OO f 

I(p) = dt d 3 r S(X>; r, t) . (Al) 

JO Jr<p 

Here the r integral runs over a sphere of radius p. Given the finite-width delta function (^j), 
we have 

3_T7 „ ,«r» „ ^ erf ( r /") 



Up) = / d 6 r V x CD x V, 

4vre Jr< P r 



1 



dSh x 



_ exlir/a) 

T> x V — y — L ^- 



47re Js 

where the integral now runs over the surface S of the sphere 
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(A2) 



Two obvious limiting cases emerge depending on the radius of the sphere. For p ^> a 
the error function inside the integral may be regarded as a constant equal to one, and the 
integral gives 

I(p) = ^-T>, p»a. (A3) 
3e 

In the contrary case /)Cawe may expand the error function as a power series in r. This 
immediately yields 

lim Up) = . (A4) 

p^O 

We have, in effect, shown that there is something akin to a delta function, albeit with a 
finite width ~ a, in the immediate proximity of the origin r = in the expression / dtS. 
Next consider an integral of the form 

1 = J™ dt J d 3 r' S(X>, r - r', t - t / )0(rr)e~ in( *~* /) , (A5) 

where the characteristic spatial and temporal scales of the function <fi satisfy 

Ar > a, At > a/c, At > fT 1 . (A6) 

Moreover, we assume that the cutoff of electromagnetic frequencies is much higher than Q, 
c/a ^> Q. With these conditions, it is possible to choose a length p such that simultaneously 
a < p <C Ar and p <^ cQ^ 1 = \/2ti. We divide the spatial integral into two regions: a 
sphere of radius p, and the complement of the sphere. In the outer region the function 5 a 
acts effectively delta function, so we have 



I r > p = — I if'r'V X[VxV) 



r't- |r-r / |/c)e ifc l r - r 'l 



47TC0 Jr>p 

f dVK(X>;r-r')0(r',t- |r-r'|/c), (A7) 

J r>p 

where K is the standard dipole radiation formula (JlQ|). On the other hand, in analogy with 
the previous expression (|A3|), the integral over the inner region gives 

J r<p = -5=-X>0(rt). (A8) 
3e 
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Summarizing, when the result is used in conjunction with a smooth function of r and t 
that also varies slowly in comparison with e~ tnt , we may write 



roo 

/ dtS(D;r,t)e~ 
Jo 

1 r Stt i 

(A9) 
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8tt 



Besides, this prescription comes with explicit directions how to handle integrals involving 
the dipole radiation. Because of the 1/r 3 divergence of the dipole radiation K, such integrals 
are generally not absolutely convergent, and their values depend on how they are done. It 
is clear from our development that the proper way to carry out such integrals is to remove a 
sphere of a finite radius p around the divergence, calculate the integral, and then let p — > 0. 
In practice, this is the same as doing the integral in spherical coordinates with the origin at 
the divergence, and integrating over the angles first. This is the prescription adopted in the 
present paper. 

One may wonder where precisely our rule for handling the 1/r 3 singularity came from, 
and whether there are plausible alternatives. Mathematically, our rule originates from the 
assumption that the photon modes were truncated in a manner that preserves the isotropy 
of photon phase space. We surmise, albeit without proof, that our prescription, indeed, is 
essentially unique if there is to be no intrinsically favored directions for photons. 

2. Vectors and matrix elements 

Alongside with the Cartesian unit vectors, we introduce the conventional circular unit 
vectors as 

e + = --^=(ei + e 2 ), e = e 3 , e_ = -^=(ei - ie 2 ) , (A10) 
These are orthonormal, in that 

e* • e CT / = <W, 5Ze ff e* = l. (All) 
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We also define a shorthand for the Clebsch-Gordan coefficients, 

(J e M; U g \laJ g m) = (M\ma) = (ma\M) . (A12) 
The dipole operator is defined as 

d = Vj2 \J e M)(M\ma)(J g m\ e* a + H.c. , (A13) 

mMcr 

where T> is the reduced dipole moment matrix element that would pertain to a transition 
with unit Clebsch-Gordan coefficient. In this paper we have chosen T> to be real. The reason 
for the complex conjugate in Eq. ( |A13|) is that we want a light field with the polarization 
e + to drive transitions with M — m = 1. The dipole matrix elements are explicitly 

d Mm = Vj2(M\ma) e* a , d mM = d* Mm . (A14) 

In particular, it may be verified from the orthonormality of the Clebsch-Gordan coefficient 
that these matrix elements satisfy 

d-Mm ■ d m M' = T^ 2 ^mm' ■ (A15) 

m 

In the special case J g = — > J e = 1 we have 

< — 1|0 - 1) = (0|00) = (+1|0 + 1) = 1 . (A16) 
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